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Abstract 

This paper discusses the spectral collocation method for numerically solving non- 
local problems: one dimensional space fractional advection-diffusion equation; and 
two dimensional linear/nonlinear space fractional advection-diffusion equation. The 
differentiation matrixes of the left and right Riemann-Liouville and Caputo fractional 
derivatives are derived for any collocation points within any given interval. The stabil- 
ities of the one dimensional semi-discrete and full-discrete schemes are theoretically 
established. Several numerical examples with different boundary conditions are com- 
puted to testify the efficiency of the numerical schemes and confirm the exponential 
convergence; the physical simulations for Levy-Feller advection-diffusion equation 
are performed; and the eigenvalue distributions of the iterative matrix for a variety 
of systems are displayed to illustrate the stabilities of the numerical schemes in more 
general cases. 

Keywords: Riemann-Liouville fractional derivative, Caputo fractional derivative, 
spectral collocation method, differentiation matrix, fractional advection-diffusion equa- 
tion. 

1 Introduction 



During the last decades, it was pointed out that the fractional calculus is more suitable to 
describe the memory and hereditary properties of materials and processes than the classical 
one, and hundreds of engineers and scientists are enrolled in the investigation of the frac- 
tional calculus and fractional differential equations due to their vast potential of applications 
in various practical fields, especially, for characterizing anomalous dynamics. In fact, these 
applications cross diverse disciplines, including physics, chemistry, biology, polymer, me- 
chanical engineering, signal processing, systems identification, control theory, finance, etc 
[16, 17, 21]. 

The space fractional derivative with order a € (1,2) can characterize superdiffusion 
of anomalous diffusion. The anomalous diffusion has the non-linear growth of the mean 
squared displacement in the course of time with the power-law form (x 2 (t)) ~ K a t a , where 
n a is the diffusion coefficient; and the superdiffusion corresponds to a £ (1, 2) being ex- 
actly the order of space fractional derivative [16]. The space fractional advection-diffusion 
equation can effectively describe superdiffusion with advection motion. This paper discusses 
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the polynomial spectral collocation method for the following space fractional advection- 
diffusion equation 



du(x, t) 
dt 



K a ^V%U m (x,t)-V a -Vu(x,t) + f(x,t), (x,t) € OX(0,T], (1.1) 



where n a > is the generalized diffusivity, u a is the velocity vector and all its compo- 
nents are positive, d = 1 and 2, respectively, corresponds to the one and two dimensional 
problems we will discuss, m is a nature number and m = 1 corresponds to the linear case, 
and u(x, i) is the unknown concentration function on the bounded domain SI = (—1, 1) , 
/(x, t) denotes a source, V is the usual gradient operator, and 

V^ = { P - 1 D^ + q Xt D?}, (1.2) 

with p, q > Q,p + q = 1, and -i-D" and Xs Df are left and right Riemann-Liouville 
fractional partial derivatives of order a (1 < a < 2) with respect to the s-th component x s 
of x, respectively, being defined by 

and x^ represents x with its s-th component replaced by £. 

The space fractional advection-diffusion equation (1.1) arises naturally in many practical 
problems, it can describe the probability distribution of the particles having advection and 
superdiffusion, such as, the nonlocal flows and non-Fickian flows in the porous media [22]. 
There are some models which are the special cases of (1.1). Chaves presents a one dimen- 
sional fractional advection-diffusion equation [3] (i.e., d = 1 in (1.1)) used to investigate the 
mechanism of the superdiffusion. A governing equation of stable random walks is developed 
in [2]. Ervin and Roop derive the steady state fractional advection dispersion equation in M. d 
via a continuous time random walk (CTRW) model [7]. For more applications of (1.1) one 
can see the recent reviews [16, 23, 17] and references therein. 

Several analytical methodologies, such as, Laplace transform, Mellin transform, Fourier 
transform, are restored to obtain the analytical solutions of the fractional equations by many 
authors [20, 16, 23, 17]. Similar to the classical differential equations, they mainly work 
for linear fractional differential equations with constant coefficients, but be invalid for the 
equations with variable coefficients and nonlinear equations; and even more the solutions 
obtained by analytical methods are usually transcendental functions or infinite series because 
of the nonlocal properties of fractional operators. So, in many cases the more reasonable 
option is to find the numerical solutions. During the recent years, much more efforts have 
been devoted to the numerical investigations of fractional integral and differential equations, 
such as, the finite difference method [5, 15, 18, 19], finite element method [4, 6, 7, 8, 24], 
and spectral Galerkin method [13, 14], etc. 

It is well known that the challenge of solving fractional differential equations essen- 
tially comes from the nonlocal properties of fractional derivatives. As a 'global' numerical 
method, spectral method seems to be a nature choice for obtaining high order numerical 
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schemes of solving fractional differential equations. Spectral collocation method has its spe- 
cial advantage for solving nonlinear problems, the aim of this paper is to design the polyno- 
mial spectral collocation algorithm for solving space fractional advection-diffusion equation 
(1.1) in one and two dimensions. 

The rest of the paper is organized as follows. In Section 2 we derive the differentia- 
tion matrixes for Riemann-Liouville and Caputo fractional derivatives of order a (n — 1 < 
a < n) within some given interval. Section 3 proposes the semi-discrete and full-discrete 
schemes for the space fractional advection-diffusion equation in one dimension and dis- 
cusses their numerical stabilities. Section 4 presents the numerical schemes for the linear 
and nonlinear fractional advection-diffusion equations in two dimensions, and the exten- 
sive numerical results are displayed in Section 5 to confirm the effectiveness of the spectral 
collocation method for solving space fractional advection-diffusion equations. 



2 Preliminaries 

This section focuses on deriving the differentiation matrixes for fractional derivatives, and 
making some remarks and presenting the propositions of the derived matrixes. This is the 
first and foremost job we need to do for applying spectral collocation method to solve frac- 
tional problems. The definitions of fractional derivatives and integral of order a(n — 1 < 
a < n) are as follows [20]: 

• left Riemann-Liouville fractional derivative: 

a D%u{ X ) = - 1 r -f- f 

x w r(n - a) dx n J a (x - 

• right Riemann-Liouville fractional derivative: 

' b 1 ' ~ T(n - a) da- J x (£ - x)<*-"+* ^ 

• left Caputo fractional derivative: 

a x w r(n - a) J a (x - £)a-n+i 

• right Caputo fractional derivative: 



• left fractional integral: 

• right fractional integral: 
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From the definitions of Riemann-Liouville derivatives, it is easily concluded that 



a D%[(x - aV] = r , r( 2+ 1} A x - ay-", 7 > -1, 
r( 7 + l-a) 

f-D 7 r(7 + i) 
M( X - b y] = l r(7 j + i 7 _ + a) J (& - xV-«, 7 > -1. 

For the left Caputo derivative of (x — a) 1 and right Caputo derivative of (x — b) 1 , we have 
the same formulae as (2.1) except for 7 = 0, . . . , n — 1, and both of them are zeros when 7 
taking these values. 

2.1 Ways to Evaluate Fractional Derivatives of Functions Approxi- 
mated by Polynomials 

There are two ways to evaluate the fractional derivatives of some function u N , being the 
polynomial approximation of u: one is to expand u by some special orthogonal polynomi- 
als first and then use the formulae of fractional derivatives of some particular orthogonal 
polynomials; another one is to do the Lagrange interpolation for u first and then derive the 
differentiation matrix based on the expanded interpolation function. 

For the first way, we can apply the following formulae for Riemann-Liouville fractional 
derivatives [1] 

(i + x) ^ Pt^M = r(* + q + i) f (1 + fl^tt) (8 _ a-id* 
1 j pi^ s - a) (i) "r(7 + i)r(a)A (i 4j p(^(i) u xj ^ 

(2.2) 

where 7, 5 > — 1, a > 0, — 1 < x < 1, and Pn' 5 \x) is the Jacobi polynomial of degree 
n with respect to the weight function w(x) = (1 — .t) 7 (1 + x) s . Applying the properties of 
Riemann-Liouville fractional derivatives a D" a D~ a = I and x D^ x D^ a = I, the above 
formulae can be rewritten as 

-x^((1 + x)^P(— '+«>(*)) = r( " + ^" + 1) (l + ,)^)(,) 



r(n 

B Uf ((1 - xr +a P^ s - a \x)) = £^±^±^±H(l - xfP^ 5 \x) 



(2.3) 



This means that the functions {(l+x) s+a P^ a ' S+a) (x)} or {(l-x) 7+a P^ 7+Q ' ,5_Q) (x)} 
can be chosen as bases to expand the function u. Actually, these basis functions are some 
form of the generalized Jacobi polynomials given in [10], which form a complete orthogonal 
system in some weighted L 2 ([— 1, 1]) space. If the fractional derivative in (1.1) is one sided 
Riemann-Liouville fractional derivative, namely, q = or p = in (1.2), these techniques 
work very well; we can choose {(1 + x) s+a pj{ y ~ a ' S+a \x)} as basis functions to expand 
u when q = and use {(1 — x) 1+a P^ y+a ' S ~ a \x)} as bases when p = 0. 
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For the general form of (1.1), including one sided (p ■ q = 0) and two sided cases 
(p ■ q ^ 0), we have to seek the other way to derive the differentiation matrix. First we do 
the Lagrange interpolation of u(x, t) at Legendre-Gauss-Lobatto/Chebyshev-Gauss-Lobatto 
points {xi : i = 0, . . . , N}, and expand it as a power function of (x — a) ( or {x — b) ); 
then using (2.1) leads to the evaluation of the fractional derivatives of u N . The Legendre- 
Gauss-Lobatto and Chebyshev-Gauss-Lobatto points are the scaled (from [—1, 1] into [a, b]) 
roots of the polynomials (1 — x 2 )P N (x) and (1 — x 2 )T N (x) respectively, where Pn(x) 
and Tjv(x) are the Legendre and Chebyshev polynomials of degree N respectively. Gener- 
ally, we need to use numerical methods to get the Legendre-Gauss-Lobatto points; but for 
Chebyshev-Gauss-Lobatto points, there exists the formula Xi = — cos(^-), i = 0, . . . , N 
in the interval [—1,1]. The Lagrange interpolation of u(x, t) is 



N N 



,{x,t) = Y J <x l ,t)k{x), h{x)= \\ - _^ - , (2.4) 

i=0 J=0,i#i 



and the Lagrange polynomials can also be written as 

[x - Xi)Q' (xi)' < ' 2 ' 5 ' 1 

where Q(x) = (1 — x 2 )P N (x) for Legendre-Gauss-Lobatto points or Q(x) = (1 — 
x 2 )T N (x) for Chebyshev-Gauss-Lobatto points. In the next subsection, we present the de- 
tails of deriving the differentiation matrixes and discuss their properties. 



2.2 Differentiation Matrixes for Fractional Derivatives 

Although we choose Legendre-Gauss-Lobatto/Chebyshev-Gauss-Lobatto points as the col- 
location points for designing the numerical schemes. The differentiation matrixes will be 
derived on any collocation points {xi : i = 0, . . . , N} within the given interval [a, b] for 
the left and right Riemann-Liouville and Caputo fractional derivatives of order a (n — 1 < 
a < n), where N > n is assumed. First we use the method given in [11] to expand k(x) as 



N N 
j=0,j=£i Xl Xj k=0 

The idea of the method is briefly described as follows. The nodal polynomial can be ex- 
panded as 

m 

p(z) = l[(z - Zi) = z m + dz"- 1 + • • • + cm. (2.7) 

i=l 

Consider the reciprocal polynomial 

q(z) = z m p{-) = TT(1 - zzi) = 1 + Cl z + ■ ■ ■ + c m z m , (2.8) 
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and 



oo oo m 



E Y^t. = - E * E ^ = - E (E ^V" 1 - cu<*> 



<z'M = E(" z ») II (i-^i) =ci + 2c 2 z + ---+mc m z m - 1 , (2.9) 

i=l j=lj^i 

then 

g'(*) 

i=1 J=0 fc=1 i=1 

m 

Denoting = X! z f > men me coefficients Cj can be determined by comparing the coeffi- 

i=i 

cients of the above equality, which obtains 

c = 1, Cj = — r(s<co + Si-ici H h sicj-i), i = 1,2, . . . ,m. 

i 

Remark 1. Due to the limits of the machine precision, the errors of the coefficients in (2.6) 
calculated by computer may explode with N being large. For N = 15, 20, 25, 30, 35, 
the max errors of the coefficients with respect to the Legendre-Gauss-Lobatto points are 
1.30 x 10~ 8 , 1.46 x 10~ 5 , 9.28 x 10~ 3 , 10.33, 5304. In order to ensure the numerical 
accuracy, we can use multiple precision packages, e.g., the Multiple Precision Toolbox for 
Matlab[25]. 

Remark 2. The coefficients in (2.6) can also be obtained from a Vandermonde system which 
is formulated by substituting the N + 1 collocation points Xi into (2.6). 

Thanks to the equalities (2.1), i.e., a D%[(x - a)f] = r ^+t-i) ( x ~ a ) 7 ~" and ( Z6 )' 
the left Riemann-Liouville fractional derivative of u N (x, t) in (2.4) becomes 

N 



aD" U N ( x , *) = aD" ^ u(xi, t)k(x) 
i=0 

N N 

= Y, u(xi, t) a D° ( Y ^,k{x - a) N ~ k ) (2.1 1) 



i=0 fe=0 
N N 



'T(N -k + l-a) 
i=0 k=0 y ' 

Then the differentiation matrix lD^ 01 ' of the left Riemann-Liouville fractional derivative is 
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For the right Riemann-Liouville derivative of u N (x, t), it follows from (2.1) that 



x x 



\N-k 



3=0,3^1 J 
N 

x D%u N {x,t) = x D%^2,u(xi,t)li( 

N 

= ^2u(xi,t) x D%(^2<l itk {x - b) 

i=0 k=0 y ' 

Thus the differentiation matrix of the right Riemann-Liouville fractional derivative 



1-0 (2.13) 

N N 

\N-k 



is 



*Dtf = E^ ( X- r fc + l ^ V ^- fc -, i,j = 0,...,N. (2.14) 

k=0 V ; 

By the same approach for Riemann-Liouville fractional derivatives and notice that for n — 

1 < a < n, 

S f DJ[(aj-o)T] = °Dg , [(a-6)*' r ]=0,7 = 0,...,n-l ) (2.15) 
we obtain the differentiation matrixes of the left and right Caputo fractional derivatives, 

g°ir' = E^ n 7^1- a) ^ -»"""•• < il6 > 



c D M ;^ d (-l)»- k T(N-k+l) N _ k _ a 
R^ji-l^^k r{N _ k+1 _ a) P *J . M-0,...,JV, 

(2.17) 



Remark 3. The norm of the differentiation matrix is infinity if some point Xj = a is 

used, and the norm of rD^"' is infinity if some point Xj = b is chosen, this is because that 
the Riemann-Liouville fractional derivatives a D" [(x — a) 7 ] and X D^ [(b — x) 7 ] in (2.1) are 
infinity at x = a and x = bfor 7 < a respectively. But this situation does not occur in the 
differentiation matrixes of the left and right Caputo fractional derivatives. 

Remark 4. For a —> n, the matrix ^r)( a ) in (2.16) is the differentiation matrix of the integer 
derivative of order n; particularly the differential matrixes (2.22) and (2.23) of the first order 
derivative on Legendre-Gauss-Lobatto and Chebyshev-Gauss-Lobatto points respectively, 
used in the numerical algorithm in the sequel, can also be obtained by (2.16) indeed. It 
coincides with the property that ^ D"f(x) = ■4-^f(x) when a —> n [20]. 

Proposition 1. If the collocation points {xi : i = 0, . . . , N } within the interval [a,b] 
satisfy Xi — a = b — x^-ifor i = 0, . . . , N, then for the coefficients in (2.6) and (2.13) 
and the differentiation matrixes of fractional derivatives, the following holds 
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(1) c^ k = {-l) N + k d N ^ k , i,k = 0, . . . ,N; 

(?) r D ( "'- D nl") C r)( a ) — C r)( a ) ■ ■ r, at. 

\4) L^ji —RU N _j N _ i , L u 3 i R N—j N — i' 'iJ U, . . . , iy , 

in fact, the second item (2) says that the differentiation matrix of the right Riemann-Liouvillef Caputo) 
fractional derivative comes by a 180 degrees rotation of the differentiation matrix of the left 
Riemann-Liouville( Caputo) fractional derivative. 

Proof. From the equality 

N N , > , s N 

it is easily observed that 

From (2.13) and Xj — a = b — xjy-j, it yields 



Zjv-i(a;) =J2 d x-i*(x-b) N - k = 



k=0 j=Q,jjtN-i XN ~ i x i 



N 

n 



(x-b)- (x N -j - b) 
j=o,&i ( XN ~ l ~ b ^~ ( XN ~i ~ b "> 

L (xi - a) - (xj - a) 



(2.20) 



AT 

Comparing the coefficients, we obtain c^k = {—^) N+k dN-i.k- And by (2.12) and (2.14), 
we have 

ji? W-V c , r(JV-fc + l) o) AT-fc-q 

* "to ' r(iV-fc + i-«) 1 J j 

-YV W-k + i) N . k . a (2.2D 

- 

— RMjy-j N-i- 

Using the similar way as above, we can also prove £-Dj" = Si^N-j N-i- ^ 
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In [12], the differentiation matrixes based on Legendre-Gauss-Lobatto and Chebyshev- 
Gauss-Lobatto points of first order derivative are, respectively, given by 

N(N+ 1) 



D L = 



and 



where 



PN{.Xi) Xj — Xi 

0, 

N(N + 1) 
4 ' 

2N 2 + 1 



* = J = 0, 
i + j, 

i = je[l,...,N-i], 
i = j = N, 
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Ci X j ' " X 'i 
X ; 



2(1 
2N 2 + 1 



G 



i = j = 0, 

i = je[l,. 
i=j = N, 



,N-1], 



(2.22) 



(2.23) 



2, i = 0,N, 
l, ie[l,...,JV-i]. 

In the following, we will also use the denotation Dji which means both D 1 ^ and work 
over there. 



3 Spectral Collocation Method for Space Fractional 
Advection-Diffusion Equation in One Dimension 

We consider the one dimensional linear case of (1.1) with the Dirichlet/Neumann/mixed 
boundary conditions, i.e., the following space fractional differential equation 

du(x,t) 



dt 



K a V a u(x, t) - v a Vu(x, t) + f(x, t), - 1 < x < 1, t > 0, 



du(— 1, t) , . 
a lU (-ht)-Pi Qx =9i(t), t>0, 

a 2 u(l,t) + p 2 =92(t), t>0, 

u(x,0) = h(x), -1<x<1, 

where 1 < a < 2, n a , v a > 0, and a 2 + 01 > 0, a*, ft > 0, i = 1, 2; ft = fa = 
corresponds to the Dirichlet boundary condition, a.\ = u-i = the Neumann one, and 
otherwise the mixed one. We will apply polynomial spectral collocation method based 
on Legendre-Gauss-Lobatto/Chebyshev-Gauss-Lobatto collocation points to approximate 
(3.1). 
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3.1 Spectral Collocation Method 

; polynomial s 
ce of polynom 



The polynomial spectral collocation method for problem (3.1) is to find u N £ Pat, being the 
space of polynomials of degree equal or less than N, such that 



= K a V a U N (Xj,t) - UaVu N (Xj,t) + f(Xj,t), 1 < j < N - 1, 

^(-l.')-A^#^=^), (3 . 2) 

a 2 u N (l,t)+(3 2 du ^ t] =g 2 (t), 
ji N (xj,0) = h(xj), 0<j<N. 
That is, for the interior points Xj, j = 1, . . . , N — 1, excluding the points —1 and +1, 

d "*ffi' f) = £ (^(p^lf + g - i/al^M^t) + /(*,-, t), (3.3) 

i=0 

and for boundary and initial conditions, the following hold 

JY 

a\u N {-l,t) - Pi^D 0i u N (xi,t) = 5i (t), 
i=0 

* (3.4) 

a 2 M JV (l,t) + P 2 2_^D Ni U N {Xi,t) = 32(0: 
i=0 



N (xj,0) = h(xj), j = 0,...,N. 



Thus, we get a 2-by-2 system for the computation of the boundary values 

N-l 

(an - f3iD 00 )u N (-l,t) - /3iD 0N u N (l,t) = gi(t) + /?i ^ D ai u N (xi,t), 



P 2 D NQ u N {-l,t) + (a 2 + (3 2 D NN )u N (l 1 t) = g 2 (t) - j3 2 ^ D Ni u N {x il t). 

i=i 

There are several ways to discretize (3.1) in the time direction, here we take the 9 scheme. 
Then the fully discrete scheme for (3.1) is to find u N G Ptv such that for 1 < j < N — 1, 

' D T U k+1 { Xj ) = Ka V a U k+1 ( X] ) - v a VU k+1 ( Xj ) + f h g + \x ), 
r)77 fe (— 1 ) 

ai U k (-l)-^ ° U „ { ' =gx{t k ), 

, 0x (3.6) 

a 2 U k (l)+(3 2 ^y=g 2 (t k ), 

U Q {xi) = h{xi), i = 0, . . . ,N. 

where U k ( X j) = u N (xj,tk), f k (xj) = f(%j,~tk), tk = kr, r is the time step size and k 
is an integer, < < 1, and the notations D T U k+1 and v$ +1 are used as 

D T U k+1 ( Xj ) = Uk+1 ( x i)- uk ( x i) ; v k+i = 6v k+x + (1 _ e y_ 
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3.2 Stability 

Taking Legendre-Gauss-Lobatto points as collocation points, we carry out the stability anal- 
ysis of (3.6) with the Dirichlet boundary condition, i.e., a,\ = = 1 and f3i = = 0, 
and naturally we can take gi(t) = gzit) = 0. Generally, the discrete inner product and 
norm are defined as follows, 

N 



{u,v) N =^u(x 1 )v{x l )w i , \\u\\ N = y (u,u) N , (3.7) 



2 = 



where {x{} and {wj are the Legendre-Gauss-Lobatto points and the corresponding quadra- 
ture weights given by 



2 



(l-xi)P N (xi) = 0, Wi = — — (Pjyfo)) *, 0<z<A, 



-2 



and —l = XQ<X\<---<Xi<---<x N = l. 

It is obvious that ||</)||l 2 = ||^|| N if is a polynomial of degree less than N . And for 
any polynomial <fr £ IV > the discrete norm is equivalent to the L 2 norm [9], namely 



<\\4>\\ N <\/2 + ±\\ct>\\ L z. (3.9) 

If N is a fixed positive integer and u N G P^, being the space of polynomials of degree equal 
or less than N and with zero boundary values, then it can be expanded as a combination of 
Legendre orthogonal polynomials {P m (x) : m = 0, . . . , N}, 

u N (x) = c P Q (x) + CiP^x) + ■ ■ ■ + c N P N (x), (3.10) 

and the left Riemann-Liouville fractional derivative of u N can be expressed as 

-iD*u N {x) = c Q ip (x) + ciipi(x) H h c N ip N (x), (3.11) 

where ip m (x) = - 1 D^P m (x). 

As {P m (x)} are Legendre orthogonal polynomials, and with (3.9), it yields 



h N \\l>\\u N \\h = \\J2c m P m (x)\\h 

m=0 

N N 

C 2 m \\P m (x)\\ 2 L2 > ^— y C "n\\Pm( 



(3.12) 



m—0 m—0 

For each m, the summation J^—i 1 ^mC^i)^ * s a fixed positive number determined by N, 
and by using (3.9) and \\P m {x)\\ L 2 = (to + 1/2)" 1 / 2 (see [9]), it yields that UPmC^L is 
positive and bounded. Then there exists a constant C dependent on N such that 

N-l 

X>™(zi)wi<C1|P m (x)||*, (3.13) 
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where {xi}fL and {u)i}fL are given by (3.8). 

Thus with the inequalities (3.12), (3.13) and (£m=o a m f < (N + 1) £™ =0 we 
have 



iV-l JV iV-l 



1=1 m=0 1=1 (3.14) 

AT 

<C(iV)^c^||P m (x)|| 2 N <C(7V) 



l«X- 



m=0 



The similar result can be obtained for the right Riemann-Liouville fractional derivative of 
u N , thus it yields 



JV-l 



£(VX(*0) 2 Wi < C(N)\\u N \\l- (3.15) 
i=i 

Based on the usual weak formulation of (3.1), the semi-discrete scheme for (3.1) with the 
homogeneous Dirichlet boundary conditions reads: Find u N <E V® N such that for any v G 
there exists 

(~m~ ,V ^ N = K *(V a u N ,v) N ~v a (Vu N ,v) N + (f,v) N , 
u N (xj,0) = h(xj), j = 0,...,N, 
where Xj, j = 0, . . . , N, are Legendre-Gauss-Lobatto points. 

Theorem 1 (Stability of the semi-discrete scheme). For any fixed N, the semi-discrete 
scheme (3.16) is stable. 

Proof. Taking v — u N in (3.16), we obtain 

7^IKII« = K*(V a u N ,u N ) N - v a (Vu N ,u N ) N + (f,u N ) N . (3.17) 
Since u N is a polynomial of degree N and u N (±1, t) = 0, it has 

{Vu N ,u N ) N = ± l (ul^,t))^=^( U 2 N (l,t)-ul(-l,t))=0. (3.18) 

As u N (±1, t) = 0, by using the Cauchy-Schwartz inequality and (3.15), we have 

(V a u N ,u N ) N <C(N)\\u N \\l. (3.19) 
The above discussion implies that the following inequality holds for every fixed N, 

^\\u N \\l<C\\u N \\l + \\f\\l, (3.20) 

then from the Gronwall inequality, we have 

CT / f T \ 1/2 

\\u N (;T)\\ N < exp(— )(K(-,0)|| w + J a ll/M)H 2 N dt) • 

□ 
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Next we consider the full-discrete scheme for (3.1) with homogeneous Dirichlet bound- 
ary conditions, for any k > and U k+1 € P% such that 

U k +\ Xj )-U k ( Xj ) = KaVau k+i {xj) „ aVU K+i {xj) + f*+\ Xj ), (3.21) 

T 

where Xj, j = 0,...,N, are Legendre-Gauss-Lobatto points. The stability of the full- 
discrete scheme is given by the following theorem. 

Theorem 2 (Stability of the full-discrete scheme). For any fixed N and i < 6 < 1, the 
full-discrete scheme (3.21) is stable, and 

\\U M \\ N <cM^)(\\UY N +rf:\\f k \\l) 1/2 , (3.22) 

fe=0 

where T = Mr. 

Proof. We can rewrite (3.21) in the following form, 

( Uk+1 ~ U \ v) N = n a {^ ,v) N - v a (VU k+ \v) N + (/*+*, v) N ■ 0.23) 

Setting v = Ug +1 in (3.23) and using the Cauchy-Schwartz inequality, for | < 6 < 1, we 
obtain 

( Uk+1 - U \ U^) N = l(9\\U k+ X + (1 - 26)(U k +\ u k ) N (1 - 6)\\U%) 

>^(\\u k+1 \\ 2 N ~\\u k \\l). 

Using (3.19), for any fixed TV, it follows 

(v a u k+ \u k+1 ) N < c\\u k+l \\l < c(e\\u k+1 \\l + (i - e)\\u k \\l). 

Since Ug +1 is a polynomial of degree N and Ug +1 (±1) = 0, then 

(vc/ e fc+1 , uj; +1 ) N = \ f w(u^(0) 2 ^ = \{u k+ \i)f - \{u k +\-i)f = o. 



For the last term, 



u k+ ^ N < 9 -\\u k ^t + l -^wu k wi + 9 -wf k+ x + 

Combining the above equations, we get 

-L(\\uk+X - m x) <c( e -\\u k +X + l -^\\u k \\l) 



l\\f k+1 \\l + ^-\\f k \\l 



(3.24) 



Lastly, summing (3.24) from k = to k = M — 1, it obtains 

M M 

\\U M \\l < \\U°\\l + Cr£ \\U k \\l + r£ \\f%- (3.25) 

k=0 k=0 

By the discrete Gronwall inequality, we obtain (3.22). □ 
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4 Space Fractional Advection-Diffusion Equation in Two 
Dimensions 



4.1 Linear Case of (1.1) 



We establish the numerical scheme for the linear space fractional advection-diffusion equa- 
tion in two dimensions, 



du 

— = Ka^x + ^ y 

u(x,y,t) = g{x,y,t) 
u(x,y,0) = h(x,y), 



^{V" + ^y)u-{v x a y a Y ■ Vu + f(x,y,t), (x,y)en,t>0, 

(x,y) edn,t> 0, (4J) 
(x,y) g 0, 



where SI = (—1, l) 2 , n a > 0, v\ — = v a > 0. The Lagrange interpolation of u at the 
mesh grid points {(xi, yj) : < i, j < N} is 



N 



u N (x,y,t) = u(x i ,y j ,t)l i (x)lj(y), 



(4.2) 



where li (x) and lj(y) are the Lagrange interpolation polynomials, and Xi and yj are Legendre- 
Gauss-Lobatto/Chebshev-Gauss-Lobatto points. 



0.5 



0< 



-0.5. 



-1 -0.5 0.5 1 

Figure 1: The Legendre-Gauss-Lobatto collocation grid for N = 10. 

The spectral collocation method to problem (4.1) is to find u N <G Pjv(fi) such that 

( du N (Xr,y s ,t) = + ^ ^ t) _ ( ^ ^ )T _ v ^ ( ^ ^ t) 

+ f(x r ,y s ,t), l<r,s<N-l, 
u N (x r ,y s ,t) = g{x r ,y s ,t), r = 0,N{0 < s < N) ands = 0,JV(0 < r < N), 
k u N (x r , y s ,0) = h(x r ,y s ), 0<r,s<N. 

(4.3) 
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Furthermore, for points {(x r , y s ) : 1 < r, s < N — 1}, 

— - ^ = 2_, ( K a(P L D ?i + 1 R D ri) ~ v a D ri )u N (xi,y s ,t) 

(4.4) 

+ 2^ (« Q (p lD% + q RDfj) - v a D sj )u N (x r ,y 3 ,t) 

+ f(x r ,y s ,t). 

The values of u N on boundary points in (4.4) are known, then it leads to an ordinary differ- 
ential equations of the (N— l) 2 unknown variables. If the unknown variables are arranged in 
the sequence as [u M , . . . , mat-i,i, . . . , wi, s , . . . , u N -i, s , . . . , Ui tN -i, u W -i,w-i] T , 
where (x r , y s ,t), then the coefficient matrix M of the right hand side is 

M = K a (p(I ® D L + D L ® I) + q(I <g>D R + D R ® I)) - v a {I &D + D&I), 

where {D L ) tJ = L D$\ {D R ) t] = R D[f , {D) t] = i, j = 1, . . . , JV - 1, and 
L D (a \ R D^ and D are given by (2.12), (2.14) and (2.22), (2.23), respectively. lis a unit 
matrix of order N — 1, and ® stands for Kronecker product. 

Using 9 scheme in the time direction for (4.3), then the full-discrete scheme of (4.1) is 
that for 1 < r, s < N — 1 the following holds 

D T U k+1 (x r ,y s ) = K a (V« + V«)U k+1 (x r ,y s ) - v a \7U k+1 {x riVs ) + f$ +1 (x r ,y s ), 
U°(x r ,y s ) = h(x r ,y s ), 

and U k {±l lVl ) = g(±l, yi ,t k ), U k {x t7 ±l) = g(xi,±l,t k ), i = 0,...,N, where 
U k {x ri y s ) = u N (x r ,y s ,t k ), f k (x r ,y s ) = f(x r , y s , t k ), r is the time step size and t k = 
kr, and the notations D T U k+1 and v g +1 are used as 

D T U k+1 (x r ,y s ) = U k+ \x r ,y s )-U k (x r ,y s ) ^ fe+1 = ^ k+l _ fc _ 

T 

4.2 Nonlinear Case of (1.1) 

Our concentration in this subsection is to apply the polynomial spectral collocation method 
to seek the numerical solution of the nonlinear advection-diffusion equation 

K, v a ) T ■ V« + f{x, y, t), (x, y) G O, t > 0, 

{x,y) edn 7 t> 0, (4 ' 5) 
(x,y) G Q, 

where fi = (—1, l) 2 and the notations are the same as those in (4.1). The spectral collo- 
cation method of (4.5) is to find u N G Pn{Q) such that for points {(x r , y s ) : 1 < r, s < 



^ = Ka(V£ + V°X 

u(x,y,t) = g{x,y,t), 
u{x,y,0) = h(x,y), 
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N — 1}, there exists 



du N (x r ,y s ,t) 
dt 

i=0 



^2 (K a (p L D% + q R Dri)u™{xi,y s ,t) - v a D ri u N (x h y 8 ,t)j 



N 

+ ( Ka (P lD % + 1 R D ™i) u ™( x r,yj,t) - v»D sj u N (x r ,yj,tyj 

3=0 

+ f(xr,y s ,t), 

(4.6) 

and the boundary and initial conditions 

u N (x r ,y s ,t) =g(x r ,y s ,t), r = 0, N (0 < s < N) or s = 0, N (0 < r < N), 
u N (x r ,y s ,0) = h{x r ,y s ), 0<r,s<N, 

where x r , y s , r, s = 0, . . . , N are the Legendre-Gauss-Lobatto/Chebyshev-Gauss-Lobatto 
points. Using the Crank-Nicholson scheme in (4.6) leads to a nonlinear system 

F(U k+1 ) ={E - Z.Dj)U k+1 - ^D F (U k+1 ) m -(E + T -Dj)U k - T -D F (U k ) m 

-^(/ fc+1 + / fe Ho, 

(4.8) 

where 

U k = [Ui tl , . . . , U N _ lf i, . . . , U ljS , . . . , U k v _ l s , . . . , Ux,iV-lJ ■ • ■ ' u N-l,N-l\ T i 

and 

fk r r k fk fk fk fk fk iT 

J — • • • 5 JjV-l,lJ • • • ) Jl,si • ■ • > JN-l,si • ■ • > Jl,JV-l) • • • i JN-1,N-11 i 

fr,s =f(Xr,y S ,tk) + [Ka(Pi^rt + Q R D r0 ) U ™ ( X 0> V*> *k) ~~ ^"-DrO (- T 0, 2/s, *)] 

+ [K a (pLD^ + q R D^)u™(x N ,y s ,t k ) - v a D rN u N (x N ,y s ,t)] 
+ [n a {p l d[q' + q R D$)u™(x r ,y ,tk) -v a D sQ u N (x r ,y ,t)] 
+ [n a {p lD^ + q RD^)u^:(x r ,y N ,t k ) - v a D sN u N (x r ,y N ,t)], 

tmdD F = K a (p(I®D L +D L ®I)+q(I®D R +D R ®I)) , Dj = -v a {I®D+D®I),wad 
{D L ) l3 = L D\f , (D R ) l3 = R D\f , (D)ij = D ijt i,j = 1, . . . , JV - 1, and E = I ® I, 
I is a unit matrix of order iV — 1, and ® stands for Kronecker product. lD^"' , R D^ and 
.D are given in (2.12), (2.14) and (2.22), (2.23), respectively. 
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5 Numerical Examples 



To test the efficiency of the spectral collocation method for the fractional advection-diffusion 
equation (1.1), several numerical examples are presented. As it is difficult to find an analytic 
solution of fractional differential equations (1.1), in the following examples we assume that 
some analytic function u satisfying the given boundary conditions is the solution of (1.1), and 
the values of the (to be determined) source function / at collocation points are calculated 
by numerical method in our implementation. The key task to compute the source / is to 
calculate the Riemann-Liouville fractional derivatives of order a (n — 1 < a < n) of a 
given function v(x), and this question can be dealt with by the following formulae 

h T (k + l-a) + T(n-a)J a (x - a ~ n+1 ^ 

n « „_^ (-l)M fc >(6)(6-x) fc -" , (-1)" f b v^(0 
X b 1 J ~ fe tu T(k + l-a) + T(n-a)J x « - xj-n+l ^ 

the second term of the right hand side of (5.1) valued at collocation points can be deter- 
mined by numerical integration, and the Gauss-Lobatto-Jacobi quadrature rule is used in our 
numerical experiments. We choose the Crank-Nicolson scheme for the discretization in the 
time direction, i.e., 9 = 0.5, in all the numerical examples given in this section. 



5.1 Examples for One Dimension 

Four numerical examples in one dimension are given in this section. In order to verify the 
accuracy of the space spectral approximation with less cost, the analytical solution of our 
problem we select is a second order polynomial with respect to the time variable t, and the 
Example 1 is in such case. The numerical results are measured in the L°° and L 2 norms 
defined by 

L °° = o<fc<JV l^ 3 *' _u ^^ fe '')l' 

/A 2 \V2 

L 2 = yi^\u{x kl -) -u N (x kr )\ u k J . 

Example 1. Taking u(x,t) = (t 2 + l)(cos7ra+l) as the solution of (3.1 ) with homogeneous 
boundary conditions since u(±l, t) = 0. The coefficients in (3.1) are chosen 
1, p = q = 0.5. 

Figure 2 shows the L°° and L 2 errors to Example 1 approximated on Legendre-Gauss- 
Lobatto points at t = 1 for a = 1.1, 1.3, 1.7, 1.9 with r = 0.1, from which we observe 
that the Legendre spectral collocation method achieves high accuracy for our problem, and 
the solution converges exponentially. 

The L°° and L 2 errors to Example 1 approximated on Legendre-Gauss-Lobatto and 
Chebyshev-Gauss-Lobatto points at t = 1 for a = 1.5 with t = 0.1 are presented in 
Table 1, from which we see that the two kinds of collocation methods have almost the same 
efficiency. Figure 3 shows the eigenvalue distribution of the iterative matrix of the full- 
discrete scheme to Example 1 for a = 1.5 with r = 0.1 and TV = 6, 12, 18, 24, the 
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Figure 2: L°° and L 2 errors to Example 1 approximated on Legendre-Gauss-Lobatto 
points att = 1 for a = 1.1, 1.3, 1.7, 1.9 with r = 0.1. 



information the figure tells us is that the method we use for (3.1) is stable in the numerical 
trial. 

Next we give an example with the mixed boundary conditions. 

Example 2. Taking u(x,t) = exp(|) sinx, n a = v a = I, p = q = 0.5, and it is a 
solution of (3.1) with the following boundary conditions. 

u(-M) - qI = fflq>(|)(cos(l) -sin(l)), t > 0, 

u(l,t) + = exp(|)(cos(l) + sin(l)), t > 0. 

ax 2 

In Figure 4, we plot the L°° and L 2 errors to Example 2 approximated on Legendre- 
Gauss-Lobatto and Chebyshev-Gauss-Lobatto points at t — lfora= 1.1, 1.3, 1.5, 1.7, 1.9 
withr = 1/10000. 
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Table 1 : The L°° and L 2 errors to Example 1 approximated on Legendre-Gauss-Lobatto 
and Chebyshev-Gauss-Lobatto points at t = 1 for a = 1.5 with r = 0.1. 



Legendre collocation method Chebyshev collocation method 



N L°° 


L 2 




L 2 


6 4.69749E-03 


4.00738E-03 


7.10790E-03 


6.74624E-03 


8 9.66278E-05 


8.76620E-05 


1.60074E-04 


1.49439E-04 


10 1.67131E-06 


1.41412E-06 


2.17903E-06 


2.34726E-06 


12 1.97558E-08 


1.71243E-08 


2.66105E-08 


2.75979E-08 


14 1.86275E-10 


1.60217E-10 


2.55300E-10 


2.51841E-10 


16 1.46372E-12 


1.19140E-12 


1.79889E-12 


1.83487E-12 


18 1.46549E-14 


9.88058E-15 


1.11022E-14 




20 4.66294E-15 


3.74558E-15 


5.77316E-15 


5.14799E-15 


Example 3. We approximate the Levy-Feller advection- diffusion 


equation [15] 




- Vm(j;, i), 


- 1 < 


x < 1, t > 0, 


u(-l,t) = u(l,t) = 0, 




t > 0, 




,tt(x + 1) 
u(x 7 0) = sm( 1 2 j 


), 


- 1 < 


a; < 1, 


w/iere V Q u = p _i-D"u + q x 


and 






sin/fa - i9)tt/2) 
sin(a7rj 


sin((a 4 
sm( 


"V 2) , W<2-o. 



In Figure 5, (a) shows the numerical solutions approximated on Legendre-Gauss-Lobatto 
points for Example 3 with r = 0.1 at t = 0, 0.2, 0.4, 0.6, 0.8, 1.0, and (b) displays the 
behavior of the solution of the advection-diffusion process for different #, From Figure 
5(b), we can see that the parameter $ plays an important role in the investigation of the 
non-Fickian transport. 



5.2 Examples for Two Dimensions 

For the two dimensional case, the errors are measured in the L°° and L 2 norms defined by 

L°° = max \u(x r ,y s , ■) - u N (x r ,y s , •) I , 

0<r,s<N 

L 2 = f 2_j \ U ( X r:Vsr) - U N (x r ,y s ,-)\ UJ r UJ s J 
r,s=0 

Example 4. Setting u(x, y, t) = (t 2 + 1) exp(a; 2 + y 2 ), n a = v a = 1, p = q = 0.5, and 
it is a solution of (4.1) with the associated Dirichlet boundary conditions. 
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(a) Legendre collocation method 



-1 -0.8 -0.6 -0.4 -0.2 0.2 0.4 
Real 

(b) Chebyshev collocation method 



Figure 3: The eigenvalue distribution of the iterative matrix of the full-discrete scheme of 
Example 1 for a = 1.5 with r = 0.1. 




Polynomial degree N 

(a) a = 1.1 



LP-L"-Norm 
LP-L 2 -Norm 
CP-L"-Norm 
CP-L 2 -Norm 




Polynomial degree N 

(b) a = 1.3 




Polynomial degree N 

(c) a = 1.5 



LP_L"-Norm 
LP-L 2 -Norm 
CP-L~-Norm 
CP-L 2 -Norm 





10 1 — 1 1 ' 1 ■ ■ — " — 1 10 1 — 1 1 ' 1 " * — * — 1 

4 6 8 10 12 14 16 4 6 8 10 12 14 16 

Polynomial degree N Polynomial degree N 

(d) a = 1.7 (e) a = 1.9 

Figure 4: L°° and L 2 errors to Example 2 approximated on Legendre-Gauss-Lobatto 
and Chebyshev-Gauss-Lobatto points at t = 1 for a = 1.1, 1.3, 1.5, 1.7, 1.9 with 
t = 1/10000 (LP: Legendre points, CP: Chebyshev points). 



In Figure 6, we plot the L°° and L 2 errors to Example 4 approximated on Legendre- 
Gauss-Lobatto and Chebyshev-Gauss-Lobatto points at t = lfora= 1.1, 1.3, 1.5, 1.7, 1.9 
with r = 0.1. Figure 7 displays the eigenvalue distribution of the iterative matrix of the full- 
discrete scheme of Example 4 for a = 1.5 with r = 0.1 and N = 10, 15, 20, it is noticed 
that the imaginary parts of the eigenvalues of the iterative matrix become smaller with a 
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t=0 




t=0.2 




t=0.4 - 


o 


t=0.6 




t=0.8 




t=1 - 



° + * 

x O + 



« i?=-0.3 

o e=-o.i5 

15=0 

+ i>=0.15 
* i?=0.3 



* 9, 



(a) a = 1.8, = 0.1, JV = 20 



(b) a = 1.6, i = 1, N = 20 



Figure 5: Approximations on Legendre-Gauss-Lobatto points for Example 3 with r 
0.1. 



being large, and the real parts move towards —1. 




10 12 14 16 18 20 



Polynomial degree N 

(a) a = 1.1 



LP-L"-Norm 
LP-L 2 -Norm 
CP-L"-Norm 
CP-L 2 -Norm 




8 10 12 14 16 18 
Polynomial degree N 

(b) a = 1.3 




Polynomial degree N 

(c) a = 1.5 



16 18 20 



LP_L"-Norm 
LP-L 2 -Norm 
CP-L~-Norm 
CP-L 2 -Norm ■ 





(d) a = 1.7 



(e) a = 1.9 



Figure 6: and L 2 errors to Example 4 approximated on Legendre-Gauss-Lobatto and 
Chebyshev-Gauss-Lobatto points at t = 1 for a = 1.1, 1.3, 1.5, 1.7, 1.9 with t = 0.1 
(LP: Legendre points, CP: Chebyshev points). 
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(a) a = 1.2 (b) a = 1.5 (c) a = 1.8 



Figure 7: The eigenvalue distribution of the iterative matrix of the full-discrete scheme of 
Example 4 approximated on Legendre-Gauss-Lobatto points for a = 1.2, 1.5, 1.8 with 

T = 0.1. 

5.3 Examples for Nonlinear case 

Example 5. Setting u(x, y, t) = (t 2 + 1) cxp(x 2 +y 2 ), K a = u a = 1, p = q = 0.5, m = 
2, and it is a solution of (4.5) with the corresponding Dirichlet boundary conditions. 

The numerical results for Example 5 are summarized in Table 2, which show that the 
spectral collocation method is also applicable for the nonlinear fractional advection-diffusion 
equations. 

6 Conclusion 

The differentiation matrix plays a crucial role in the spectral collocation method, especially 
for fractional differential equations. We have derived the differentiation matrixes of the left 
and right Riemann-Liouville and Caputo fractional derivatives and successfully applied the 
spectral collocation method to handle the numerical solution of the fractional advection- 
diffusion equation. The stabilities of the semi-discrete and full-discrete scheme in one di- 
mension are theoretically established for the Legendre spectral collocation method. The 
eigenvalue distributions of the iterative matrix for a variety of systems are computed for 
further confirming the stability of the numerical schemes in more general cases. The numer- 
ical results have shown the efficiency of the proposed methods, moreover, it can be noted 
that the numerical solution converges exponentially when the exact solution is analytic. The 
performance of the spectral collocation method for nonlinear fractional advection-diffusion 
equation is also exhibited. By simulations the interesting physical phenomena for Levy- 
Feller advection-diffusion equation are discovered. 
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Table 2: The L°° and L 2 errors to Example 5 collocated on Legendre-Gauss-Lobatto 
points at t = 1 for different a with t = 0.1. 





TV 


L°° 


L 2 


TV 




L 2 


1.1 


6 


1.446E-02 


1.219E-02 


14 


2.707E-07 


2.415E-07 




8 


1.249E-03 


1.012E-03 


16 


1.309E-08 


1.195E-08 






o.Uuzi>Uj 


/.I jog-ID 


1 Q 
1 O 








12 


5.023E-06 


4.400E-06 


20 


4.800E-09 


4.962E-09 




A/ 




L 2 


TV 




L 2 


1.3 


6 


2.135E-02 


2.027E-02 


14 


3.403E-07 


3.590E-07 




8 


1.554E-03 


1.629E-03 


16 


1.759E-08 


1.829E-08 




LU 


i i nAin ri/i 
1. lU0Ji-U4 


i 1 1 nc n/i 


1 c 




A ^ 1 In no 




12 


6.298E-06 


6.665E-06 


20 


3.232E-09 


5.049E-09 




TV 




L 2 


iV 




L 2 


1.5 


6 


2.494E-02 


2.469E-02 


16 


1.962E-08 


2.043E-08 




8 


1.961E-03 


1.984E-03 


18 


8.622E-10 


9.002E-10 




i n 

1 W 








^ AQATh 1 1 

J. L t7\JLL- 1 1 






12 


7.697E-06 


7.947E-06 


22 


1.284E-12 


1.360E-12 




14 


4.087E-07 


4.234E-07 


24 


7.194E-14 


5.820E-14 


Q 


TV 


L°° 


L 2 


TV 


L°° 


L 2 


1.7 


6 


2.689E-02 


2.790E-02 


14 


4.536E-07 


4.645E-07 




8 


2.203E-03 


2.231E-03 


16 


2.200E-08 


2.225E-08 




10 


1.471E-04 


1.492E-04 


18 


8.563E-10 


9.453E-10 




12 


8.725E-06 


8.785E-06 


20 


1.263E-10 


1.260E-10 


Q 


TV 




L 2 


TV 


L°° 


L 2 


1.9 


6 


2.523E-02 


3.064E-02 


14 


4.357E-07 


4.812E-07 




8 


2.103E-03 


2.354E-03 


16 


2.255E-08 


2.343E-08 




10 


1.436E-04 


1.559E-04 


18 


2.996E-09 


2.889E-09 




12 


8.405E-06 


9.129E-06 


20 


1.208E-09 


1.297E-09 
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